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Over the past few years multidimensional 
self-consistent plasma simulations including 
complex chemistry have been developed 
which are promising tools for furthering 
our understanding of reactive gas plasmas 
and for reactor design and optimization. 
These simulations must be benchmarked 
against experimental data obtained in well- 
characterized systems such as the Gaseous 
Electronics Conference (GEC) reference 
cell. Two-dimensional simulations relevant 
to the GEC Cell are reviewed in this paper 
with emphasis on fluid simulations. Impor- 



tant features observed experimentally, such 
as off-axis maxima in the charge density 
and hot spots of metastable species density 
near the electrode edges in capacitively- 
coupled GEC cells, have been captured by 
these simulations. 
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1. Introduction 



Low pressure (0.13 Pa to 1333 Pa), cold (gas temper- 
ature 300 K to 500 K), weakly ionized (degree of ioniza- 
tion 10''' to 10'') glow discharge plasmas are used exten- 
sively in the processing of electronic materials, 
especially for etching and deposition of thin films [1]. 
They also find application in surface modification (e.g., 
hardening, corrosion resistance) and lighting. In reactive 
gas plasmas, electrons decompose the flowing feed- 
stock gas into radicals and ions. In plasma deposition, 
radicals adsorb on the wafer surface where they react to 
deposit a thin film. The film microstructure and proper- 
ties (e.g., density, resistivity) can be influenced by 
energetic ion bombardment which occurs naturally on 
all surfaces exposed to the plasma. In plasma etching, 
radicals adsorb and react on the wafer to form volatile 
products which desorb and are pumped away by a 
vacuum system. The surface chemistry can be strongly 
modified by energetic ion bombardment. Ions bombard 
the wafer preferentially along the vertical direction. 



' Author to whom correspondence should be addressed. 



enhancing the reaction rate and inducing anisotropy 
which is critical for delineating sub-half micron patterns 
in advanced microelectronic device manufacturing. 
Controlling the flux, energy distribution and angular 
distribution of ions and neutrals bombarding the wafer is 
of paramount importance in plasma systems. Also, the 
uniformity of these fluxes over large diameter 
(> 200 mm) wafers is critical for the success of indus- 
trial plasma processing equipment. 

Unfortunately, glow discharge plasmas are extremely 
complex systems in which a plethora of interdependent 
parameters can influence the process, often in a subtle 
way. Thus, slight variations in reactor design and operat- 
ing conditions can result in large variations in system 
performance. The Gaseous Electronics Conference 
(GEC) Reference Cell was conceived to serve as a com- 
mon platform for experimental and modeling studies in 
different laboratories [2] . The Reference Cell is thought 
to be a well-characterized system in which fundamental 
studies of plasma behavior can be conducted. Experi- 
mental data obtained in the Cell are also useful for 
benchmarking plasma simulations, which in turn 
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provide insight into the plasma dynamics. This syner- 
gistic experimental-modeling approach is crucial for 
furthering our understanding of plasma systems and for 
the development of predictive simulation tools which are 
useful for the design and optimization of new reactors. 
The Reference Cell can be operated in two configura- 
tions: capacitively- or inductively-coupled plasma. Most 
studies related to the Cell so far are for the capacitively- 
coupled configuration, a schematic of which is shown in 
Fig. 1 a. The bottom electrode is usually powered with a 
13.56 MHz power supply, and is separated from the 
grounded walls by thin insulators. This system has 
azimuthal symmetry, but the area of the powered elec- 
trode is much smaller than that of the grounded 
electrode. Since the time-average sheath voltage in 
capacitively-coupled systems scales with the inverse 
ratio of the electrode areas [1], a large voltage drop 
appears in the sheath over the powered electrode with a 
much smaller voltage across the sheath over the 
grounded electrode. Sometimes, a push-pull rf drive 
configuration is used in which both the bottom and the 
top electrodes are powered with a rf voltage of the same 
amplitude but 180° out of phase. This results in a dis- 
charge that is symmetric not only azimuthally but also 
axially. The time-average sheath voltages over the top 
and bottom electrodes will then be identical in magni- 
tude. The capacitively coupled Cell, henceforth to be 



designated as the Gaseous Electronics Conference Ca- 
pacitively Coupled Plasma (GEC-CCP) Cell, is usually 
operated as a relatively high pressure (6.67 Pa to 133 Pa) 
low density (charge density < 10'' m"^) plasma (LDP). 

Recently, the Cell has been operated with an inductive 
coil to generate a low-pressure (< 6.67 Pa) high charge 
density (> 10'' m"^) plasma (Fig. lb). High density 
plasma (HDP) sources are becoming increasingly im- 
portant in microelectronics as device dimensions con- 
tinue to shrink. Low pressure provides more directional 
ion bombardment and better plasma uniformity over 
large diameter wafers. High plasma density ensures that 
the etch or deposition rate are comparable to those found 
in high pressure LDP systems. The inductively coupled 
Cell will henceforth be designated as GEC-ICP 

Most experimental studies with the GEC cell have 
used noble gas plasmas (Ar, He), because they are sim- 
pler than reactive gases (CI2, SFg). This also facilitates 
comparison of experimental data with simulation results 
since the electron impact collision cross sections are 
known and the noble gas plasma reactions are better 
characterized as compared to reactive gases. Further- 
more, noble gas plasmas are electropositive (i.e., the 
negative ion to electron density ratio is much less than 
unity, {nJne«\) and as such are simpler to experiment 
with, both in the laboratory and on the computer as 
compared to electronegative {nJne>\) plasmas. 




Fig. 1(a). Schematic diagram of the capacitively-coupled GEC Reference Cell. 
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Fig. lb. Schematic diagram of the inductively-coupled GEC Refer- 
ence Cell (Ref [68]). 



2. Problem Statement 

Given the reactor configuration shown in Fig. la or 
Fig. lb and a set of operating parameters (feedstock gas 
composition and flow rate, gas pressure, voltage or cur- 
rent waveform applied to the electrode), determine the 
following: The space and time variation of electron, ion, 
and neutral species densities and velocities, the flux and 
energy distribution of ions and neutrals bombarding the 
electrodes and their uniformity across the electrodes, 
the power deposited into the plasma, and the potential 
and current distribution in the system. If the reactor is 
loaded with a wafer to etch, one is in addition interested 
in the etch rate, uniformity, anisotropy (shape of micro- 
scopic features etched into the wafer), selectivity, and 
radiation damage. The level of detail one can obtain 
depends on the type of simulation used. For example, 
fluid simulations can't provide the species distribution 
functions but only averages over the distribution. 



3. Plasma Simulation 

Modeling and simulation of glow discharge plasmas 
has emerged as a tool for enhancing one's intuition 
about the physicochemical processes occurring in the 
plasma, for understanding the complex spatiotemporal 
plasma dynamics, and for assisting in the design of new 



reactors or the optimization of existing ones. Two- 
dimensional plasma reactor simulations have been 
reported in the literature in recent years [3-5] . However, 
these works focused on the transport and reaction of 
neutrals only (neutral transport and reaction models). 
The electron density was assumed to have a uniform or 
Bessel function profile, and the electron temperature 
was not calculated as a function of space and time in the 
reactor. These studies did not solve the problem of neu- 
tral radical transport and reaction in a self-consistent 
manner. The radical source terms (by electron-impact 
dissociation, for example) were estimated and the con- 
servation equations for mass, momentum, and energy 
transport of the neutrals were solved to obtain the fluid 
velocity profiles, neutral gas temperature and the con- 
centration distribution of radicals. Charged particle 
transport was not considered, and the effect of plasma 
gas composition (different from the feedstock gas com- 
position) on the plasma properties was not accounted 
for. 

Up until a few years ago, simulations that solved for 
the rf plasma dynamics (using the so-called glow dis- 
charge models) were confined to one spatial dimension 
(1-D) [6-19]. In addition, most of these simulations did 
not solve for the transport and reaction of neutrals. Self- 
consistent rf plasma simulations which solve for the 
coupled effects of charged and neutral species transport 
have only been reported within the past few years in 1-D 
[16,18,19] and 2-D [20-22]. Two dimensional simula- 
tions are particularly useful since they can address the 
important issue of plasma uniformity and the spatiotem- 
poral plasma dynamics along both the radial and axial 
direction. As of this writing, most 2-D simulations do 
not include neutral transport and chemistry and have 
considered noble gases (argon and helium) only, not 
reactive gas plasmas [23-27]. We are aware of only a 
few 2-D plasma simulations which couple the neutral 
transport and chemistry with the glow discharge in a 
self-consistent manner [20-22]. Several more groups are 
now working on this problem. In view of the above 
discussion, multidimensional self-consistent plasma re- 
actor simulation is still at an early stage of development. 

There are three kinds of glow discharge simulations: 
fluid, kinetic and hybrid. Fluid simulations integrate 
moments of the Boltzmann equation (see below) 
describing species density, momentum and energy con- 
servation. They require some assumptions regarding the 
species distribution function to achieve closure of the 
equations. Kinetic simulations, such as Particle-In-Cell 
with Monte Carlo Collisions (PIC-MCC) [28-30] or 
Direct Simulation Monte Carlo (DSMC) [31] follow the 
motion of a number of superparticles accounting for 
their collisions using Monte Carlo techniques. Kinetic 
simulations yield the particle distribution functions as an 
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output of the simulation. They are considered more 
accurate than fluid simulations in the sense that they can 
be used at low pressures (mean free path comparable to 
or longer than the characteristic length scale, or 
Knudsen number Kn = \IL>\) or for highly nonequi- 
librium situations. However, there is some evidence that 
fluid simulations can perform well even at low pressures 
for which they should be considered suspect [32]. 
Kinetic simulations are computationally intensive as 
compared to fluid simulations. Hybrid simulations have 
been developed [16, 22] in an attempt to preserve the 
accuracy of kinetic simulations and at the same time 
reduce the computational burden. 



4. Model Formulation 

This paper emphasizes the fluid simulation approach 
since all studies of the GEC Cell are based fully or 
largely on this approach. The reader is referred to the 
literature for examples of kinetic plasma simulations 
using PIC/MCC [28-30] or DSMC [31]. 



The right hand side of Eq. (1) is an integral containing 
/which is the unknown [33]. Eq. (1) is extremely diffi- 
cult to solve because it is an integrodifferential equation 
in three space dimensions (x, y, z), three velocity di- 
mensions {U;a Uy, u^) and time (?). Kinetic simulations 
such as DSMC and PIC/MCC solve for the species dis- 
tribution function by integrating the Boltzmann equa- 
tion using pseudoparticles. On the other hand, fluid 
simulations solve equations for "average" quantities such 
as the directed velocity. The fluid equations are obtained 
from the Boltzmann Eq. (1) by first multiplying that 
equation by (^, where (^ is a constant or a function of the 
velocity u , and then integrating over the velocity space 
to obtain the "average" of (^. The resulting equations are 
called moments of the Boltzmann equation [34]. 

When (^ = 1 , the species ( number density continuity 
equation (where ( can be electrons i = e, positive ions 
( = +, negative ions ( = -, or neutrals ( = n) is obtained 
as the zeroth moment. 



drii 
~dt 



+ V- {n,Yi) = 2 R„ 



(4) 



4.1 The Fluid Equations 

In general, species (electrons, ions, neutrals) trans- 
port can be described by the Boltzmann equation 



dt 



+ u -Vrf+h • V„/ = 



8?' 



(1) 



which is a continuity equation in phase space (r, u), 
where r is the spatial location vector and u is the species 
velocity vector. The acceleration is h = F /m , where F is 
the force acting on the species and m is the mass. For a 
particle with charge q, F = q{E + uXB), is the Lorentz 
force, where E and B are the electric field and magnetic 
induction, respectively. The right hand side of Eq. (1) is 
the so-called collision integral, which accounts for the 
effect of collisions on the velocity distribution function 
(VDF) /. The VDF is defined such that the number of 
particles diV with velocities between u and « + du in a 
given volume dr of configuration space is given by 



dA^=/drdH. 



(2) 



Consequently, the number density as a function of posi- 
tion can be obtained as 



nir) = Jfdu, 



where the integral is over the velocity space. 



(3) 



where n, and Vi are the density and directed velocity of 
species (', respectively. Ry is the rate of production or 
loss of species / due to the volumetric reaction j . 

The first moment is the (vector) momentum equation 
which can be derived by setting (p = mu , 



dt 



iniiniVi) + V • ininiiYiYi) 



-VPi + niiriihi - riiniiYiVmi, 



(5) 



where Vm, is an effective momentum exchange frequency 
of species ( and Pi is the partial pressure given by 
Pi = riiKTi, where k is the Boltzmann constant and T, is 
the species temperature. Equation (5) assumes that the 
pressure tensor is isotropic which appears to be a good 
approximation in the absence of strong magnetic fields. 
In Eq. (5), the terms represent (in order from left to 
right) local acceleration, convective acceleration (iner- 
tia), motion due to pressure gradient, motion due to the 
presence of force fields and momentum exchange due to 
collisions with the background species. 

When the terms on the left hand side of Eq. (5) are 
negligible (see Ref. [17] for a discussion of this point), 
one obtains 



Yi = Ui-E - - V {rii KTi) 

OTiV„; \ rii 



(6) 
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Since the flux for each carrier is /, = n,T, one can write 



/, =^,yA,M,£' -DiVui, 



(7) 



where /ij and D, are the particle's mobility and diffusiv- 
ity respectively, and the species temperature has been 
taken outside the differential as if the temperature were 
constant. Equation (7) is known as the drift-diffusion 
approximation and is often used in place of the full 
momentum Eq. (5) for simplicity [6-10]. The drift-dif- 
fusion approximation has been questioned at low pres- 
sures and in the sheath where species inertia may not be 
negligible [34]. Comparing Eqs. (6) and (7) one obtains 



M.- 



_2l 



m,v„ 



(8) 



Here Kj is the thermal conductivity and Hjj is the energy 
loss due to collision of typej . In order from left to right, 
the terms in Eq. (12) represent the time rate of change 
of total energy, convective transport of energy, rate of 
work done by pressure forces, energy exchange with the 
force field, energy transport by conduction, and energy 
loss due to collisions. 

An equation for the thermal energy [/, (For a 
Maxwellian distribution [/, = 3/2kT',) can be obtained 
by taking the dot product of v, with the momentum Eq. 
(5), and subtracting the resulting equation from the total 
energy Eq. (12) [34]. However, a simplification is 
frequently made in plasma simulations, namely that 
Ui »m,v?/2, and 5^ in Eq. (12) is simply replaced by 
[/,. The resulting equation, written here for electrons 
(i = e), reads. 



£>, 






(9) 



|(f«.Kr.) + v-^. 



Combination of Eqs. (8) and (9) results in Einstein's 
relation. 






(10) 



Although electrons are mobile enough to respond to the 
variations of the electric field at 13.56 MHz (i.e., 
Vme»u), and the electron inertia can be neglected), ions 
are massive and can't follow the field faithfully. Recog- 
nizing this fact, Richards et al. [8] introduced an equa- 
tion for an "effective" field to which the ions respond. 



+ eJ,-E + 3v„, ^ nMTe -T,) + 2 R^jHej = (13) 



with the total electron energy flux given by 



q, = -K,^T, + ^KTJ, 



(14) 



The thermal conductivity of a monatomic "gas" is given 
by 



dEf 
dt 



= v^i{E-Ef) 



(11) 



Ki = - KDifii. 



(15) 



That way, the drift-diffusion Eq. (7) can be used for both 
electrons and ions, except that the actual electric field is 
replaced by the effective field, in the case of ions. Of 
course, when the full ion momentum Eq. (5) is used, 
Eqs. (6) and (7) become immaterial. 

The second moment of the Boltzmann equation is the 
(scalar) energy conserwtion equation which can be 
derived by setting 4> = mu ■ u/2 = ^, where 5^ is the 
total (kinetic plus thermal) energy. For particle i, 
^ = m, vj 12 + Ui with [/, the thermal energy 



dt 



+ V • {ntViS:,) = -V • P,Y, 



+ n,mjii • V, + V • ^, Vr, - 2 RijH^j . (12) 



In Eq. (13), the electron energy loss term has been 
decomposed into losses due to elastic collisions, and 
inelastic collisions (last two terms on left hand side). 
Some authors write Eq. (13) in terms of the thermal 
energy [/,, i.e., they don't make the substitution 
[/, = 3/2kT', [27, 35]. This way, they have an equation in 
terms of the mean electron energy instead of one in 
terms of electron temperature. The two formalisms are 
equivalent. 

Depending on the approximations made, different 
sets of the equations shown above are used by different 
authors. Most often, the drift diffusion approximation is 
made for both electrons and ions [6-10, 20]. Other au- 
thors solve the full momentum equations for either elec- 
trons or ions (using drift-diffusion for the other species) 
[25,26,36] or for both electrons and ions [11]. 



477 



Volume 100, Number 4, July- August 1995 

Journal of Research of the National Institute of Standards and Technology 



4.2 Electron-Impact Rate Coefficients and Trans- 
port Properties 



kj' 



1 



aj{s)u,{s)f{s)ds. 



(17) 



The above derivation of the fluid equations tacitly 
assumes that the velocity distribution function is 
Maxwellian (introducing a "temperature" automatically 
implies a Maxwellian distribution). Although the VDF 
of ions and neutrals can be non-Maxwellian at low pres- 
sures, it is the electron energy distribution function 
(EEDF) which is of primary concern. Experimental 
measurements in capacitively coupled rf plasmas have 
shown non-Maxwellian EEDFs [37]. It is thought that 
the most important effect of a non-Maxwellian distribu- 
tion would be on the electron-impact reaction rate coef- 
ficients, especially those for high threshold energy pro- 
cesses (e.g., excitation, ionization) which are sensitive to 
the shape of the tail of the EEDF, or processes with 
sharp resonances. 

Rate coefficients for electron impact reactions are 
needed as input to the fluid equations. For example, 
when applied to electrons in an argon discharge, the 
right hand side of Eq. (4) reads. 



E 



Rej = kiNlle + killAT* rie + k^pllX,, 



(16) 



where A^, n^ and Kat* are ground state argon, electron, 
and metastable species density, respectively. The rele- 
vant reactions producing electrons and the correspond- 
ing rate coefficients are shown in Table 1 (reactions 2, 3, 
and 6). Electrons are produced by direct ionization of 
ground state Ar, by ionization of metastable atoms (two- 
step ionization) and by metastable-metastable collisions. 
The electron-impact rate coefficients are calculated by 
an expression of the form. 



where kj is the rate coefficient for process j, ctj is the 
corresponding collision cross section, /(e) is the EEDF, 
and s is the electron energy s = mjill2, u^ being the 
magnitude of the electron velocity, m^ = ImJ. In addition, 
electron transport properties (momentum exchange col- 
lision frequency, mobility, diffusivity, etc. ) are needed. 
These can also be obtained as the appropriate integrals 
over the distribution function [38]. 

Kinetic simulations can predict the EEDF and Eq. 
(17) can be integrated to calculate the rate coefficients. 
In fluid simulations, the electron-impact rate coeffi- 
cients are expressed as a function of the mean electron 
energy. Two approaches are used: (a) the EEDF is as- 
sumed, for example Maxwellian or Druyresteyn distri- 
bution and, knowing the cross sections, Eq. (17) is inte- 
grated directly to find kj as a function of average energy, 
or (b) A spatially-independent (0-D) Boltzmann equa- 
tion solver is used [either finite difference or finite 
element or Monte Carlo solution of Eq. (1)] to calculate 
the distribution function (for a given gas composition) 
and in turn the electron-impact reaction rate coefficients 
and transport properties as a function of E/N . At the 
same time the mean electron energy is calculated as a 
function of E/N. The results are combined to eliminate 
E/N and express the rate coefficients and the transport 
properties as a function of mean electron energy (or 
temperature; in the case of non-Maxwellian distribu- 
tion, the "temperature" is defined as 2/3 of the average 
energy, Te = 2s/3k). Since the latter is one of the depen- 
dent variables that is obtained as part of the solution 
[(see Eq. (13))], this representation is very convenient. It 
should be noted that this approach tacitly assumes that 



Table 1. Important collision processes in the agron discharge 



No. 


Process 


Reaction 




Hj(eV) 


Rate coefficient" 


1 


Ground state excitation 


Ar -1- e— > Ar* -1- e 




11.56 


^cx 


2 


Ground state ionization 


Ar -1- e ^ Ar* -1- 2e 




15.7 


h 


3 


Step- Wise ionization 


Ar* + e ^ Ar + 2e 




4.14 


kn 


4 


Superelastic collisions 


Ar* -1- e ^ Ar -1- e 




-11.56 


fcc 


5 


Quenching to resonant 


Ar* + e ^ Ar' + e 






i, = 2 X 10-' 


6 


Metastable pooling 


Ar* -1- Ar* -^ Ar* + 


Ar -1- e 




*:„p = 6.2 X 10-'° 


7 


Two-Body quenching 


Ar* + Ar ^ 2 Ar 






fe, = 3 X 10-'= 


8 


Three-Body quenching 


Ar* -1- 2 Ar ^ Ar2 + 


Ar 




fe, = 1.1 X 10-" 



° Rate coefficients for processes \-A are given as a function of electron energy in Fig. 2 of Ref. [19]. Units are 
cmVs except for fe, which is in cmVs. 
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the EEDF in the rf discharge of interest is similar to that 
which would be obtained in a DC swarm experiment. 
This approach apparently works fairly well [39]. A more 
complicated approach in which the time dependence of 
the EEDF is also taken into account has also been pro- 
posed [16]. The charged species mobility is usually 
assumed constant, although it can be expressed as a 
function of the mean electron energy. Finally, the 
charged species diffusivity is commonly obtained from 
the mobility and the species temperature using the 
Einstein relation, Eq. (10). 

In the hybrid approach of Kushner and coworkers 
[18, 22], the electron energy fluid equation [Eq. (13)] is 
dropped. Instead, the time-average electron-impact rate 
coefficients and transport properties are obtained by a 
2-D Monte Carlo simulation. Therefore, no assumptions 
need be made regarding the distribution function. 

Some 1-D RF [40, 41] and 2-D RF [24] and DC [42] 
glow discharge simulations applied the so-called local 
field approximation. In this approximation, the electron 
energy Eq. (13) is dropped. The electron impact rate 
coefficients and transport properties are expressed as a 
function of EIN . It is assumed that these quantities at a 
given point in the discharge and at a given time during 
the rf cycle are equal to the DC values that would be 
obtained using the same EIN as existed at that point in 
the discharge and that particular time in the rf cycle. 
This approximation does not allow for nonlocal behavior 
of the EVDF and is particularly bad at low pressure, for 
beam electrons (emitted by ion bombardment of the 
electrodes), and for describing electron transport in the 
sheath. 

Note that including Eq. (13) allows for nonlocal elec- 
tron transport to be captured since the rate coefficients 
are a function of mean electron energy and not EIN. 
And it is well known that fluid models predict the elec- 
tron energy to peak at the plasma-sheath interface, 
whereas the electric field peaks right at the electrode 
surface. Comparisons of fluid with kinetic simulations 
[33] show that fluid models underpredict the nonlocal 
electron behavior. For example fluid simulations predict 
a sharper electron energy peak compared to kinetic sim- 
ulations. 

In addition to electron transport and reaction coeffi- 
cients, one also needs rate coefficients for ion-neutral 
and neutral-neutral reactions and the transport proper- 
ties (mobility, diffusivity, etc.) of the heavy (ions, neu- 
trals) species. In fluid models, the ion energy distribu- 
tion function is usually assumed to be Maxwellian with 
a temperature equal to the gas temperature. When the 
full momentum equation (Eq. 5) for ions is solved, a 
drifting Maxwellian distribution is assumed. An ion en- 
ergy (Eq. 12) or a corresponding temperature equation 
have not been incorporated in the fluid models so far. 



but work towards that goal is in progress [43]. 

4.3 Maxwell's Equations 

Self-consistent plasma simulations must solve for the 
electromagnetic fields which develop in the reactor. 
Maxwell's equations relate the electromagnetic field in- 
tensity and flux density vectors to each other and to the 
sources of the field, the electric charges and currents. In 
differential form. Maxwell's equations are [44]: 



— -f. 



VX« = / + f 

at 



V -0 = p. 



V ■ B = 0. 



(18) 
(19) 
(20) 
(21) 



In the above equations E, D, H, B, J, and p are the 
electric field intensity, electric field flux, magnetic field 
intensity, magnetic field flux, current density, and 
charge density, respectively. These equations are not 
independent. For example Eq. (20) can be obtained 
by taking the divergence of Eq. (19) and making use 
of the equation of continuity of electrical charge, 
V*/ + dpidt = 0. The above equations are augmented 
with the constitutive relations B = fxH and D = sE 
where jx is the permeability and s is the permittivity of 
the medium. 

In the absence of any time varying magnetic fields (as 
is the case of typical GEC-CCP cells), dBldt = and 
Faraday's law Eq. (18) suggests that the electric field 
can be written as the divergence of a scalar. Thus, 



E = -vy. 



(22) 



where V is the electric potential. Substituting Eq. (22) 
into Eq. (20) and assuming a dielectric constant inde- 
pendent of position, yields an equation relating the 
Laplacian of the potential to the charge density, referred 
to as Poisson's equation. 



V'y=-p/e. 



(23) 



When solving the Poisson equation over a domain that 
includes different material properties, the dielectric 
constant should be kept inside the differential. 
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The equation for the propagation of electromagnetic 
fields is obtained from the Maxwell equations as 



V"£:- 



-/Llfi 






(jua 



dE 

dt 



(24) 



where use has also been made of Ohm's law, / = aE {a 
is the conductivity of the medium). Equation (24) is the 
three-dimensional wave equation for an absorbing 
medium. Similar equations can be derived for the other 
field quantities, D,B , and H [44]. Equation (24) can be 
solved for the fields as a boundary value problem. Alter- 
natively, the magnetic vector potential A formulation 
may be used. Based on the general notion that the diver- 
gence of the curl of a vector is zero (V- VX v = 0) one 
may write the magnetic field of Eq. (21) as the curl of 
a vector A , 

B = VXA. (25) 

Using Eq. (25), Faraday's law Eq. (18) can be written as 

Vx(e + ^]=0. (26) 



Then by using the identity VxW = 0, where V is the 
electric potential one obtains 



representation A^, = Ae""', fi^ = £e""' and /^ = Je""' 
respectively, where A, E and / are the corresponding 
(complex) amplitudes which depend only on space, and 
CO is the applied frequency. Equation (29) then becomes 



-j(oA 



(30) 



Recognizing that A^ is only a function of (r , z) Eq. (28) 
is simplified to 



r dr 



dA\ d^A 



dr 



+ ^TT + {(o as -r 

az 



'■)A = -ijJ, (31) 



where the complex current density is given by / = aE 
and a is the complex conductivity of the medium. Mak- 
ing the cold plasma approximation, the complex plasma 
conductivity (Xp is given by 



n^q 



»le(Veff + jO)) 



(32) 



where m^ and Veff are the electron mass and effective 
momentum-transfer collision frequency, respectively. 
The plasma dielectric constant s^ is obtained from 



E + ^ = -VV, 
at 



(27) 



where the negative sign is introduced on the right-hand 
side so that V becomes the electric potential in a static 
situation, when A is independent of time. 
One can now derive an equation for A , 



VU- 



(jbe 



dt' 



-fJiJ 



(28) 



where / is the current density giving rise to the electro- 
magnetic fields. For azimuthally symmetric systems for 
which the applied current has an azimuthal component 
only, Ar and A^ are zero, and only the A^ component of 
the magnetic vector potential survives. In this case one 
also has dV/d^ = and, using Eq. (27), the azimuthal 
component of the electric field is given by 



Es ■ 



dA» 
dt 



(29) 



and Er and E^ are given by the respective radial and axial 
gradients of the electric potential V. 

It is quite common to replace the time-harmonic elec- 
tromagnetic quantities A^, E^ and J^ with their phasor 



«P = So-J- 



(33) 



where So is the permittivity of vacuum. 

Equation (31) can be solved in the whole domain 
(including plasma, reactor walls, etc. ) by a finite differ- 
ence [22, 45] or finite element method [46]. Once A has 
been determined the azimuthal electric field is retrieved 
by using the following equation 



E^ = Re i-joAd'"). 



The power deposition Wis equal to 



W(r, z) = I iRelaE' 



(34) 



(35) 



The inductive power deposition given by Eq. (35) should 
be added to the right hand side of Eq. (13) in the GEC/ 
ICP reactor case. The assumption of Ohm's law implies 
that the power deposition is by collisional processes. At 
low operating pressures (< 1.33 Pa) noncollisional 
power deposition can take place [47, 48]. However, the 
same formulation as for the collisional case may be 
used, except that an effective electron collision fre- 
quency must be used in Eq. (32) [49]. 
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4.4 Boundary and Initial Conditions 

Boundary conditions are difficult to specify in fluid 
simulations, partly because the physics of the problem is 
not well understood (in which case even kinetic simula- 
tions have a problem) and partly because the resulting 
expressions may not be convenient to use for solving the 
fluid equations. Boundary conditions are needed for the 
species density, velocity, and temperature (or energy). In 
addition, boundary conditions for solving Maxwell's 
equations are necessary. 

The boundary condition on electron density takes the 
form of essentially a "mass balance" at the wall. The 
electron flux at the wall equals the thermal flux (assum- 
ing that electrons striking the wall are fully absorbed) 
minus the secondary electron emission flux. The latter 
equals the positive ion flux striking the wall /+ times the 
secondary electron emission coefficient y+. The thermal 
flux of electrons is given by ViuJA where Vt is the ther- 
mal velocity, equal to {?>kTJ Trni^Y^ , and «e is the elec- 
tron density, both evaluated at the wall. Therefore, 



TTTOe 



Me«-r+/+ 



(36) 



Combining Eq. (36) with the expression for the elec- 
tron flux Je as given by Eq. (7), one obtains a rather 
complex expression in terms of electron density and 
temperature, ion density and the species mobility at the 
wall. Some authors have used the simpler condition 
Me = 0. However, both kinetic and fluid simulations have 
shown that a substantial electron density can exist near 
the electrode, especially during the anodic part of the rf 
cycle [33]. Eq. (36) is written here for one type of ion 
striking the wall. For a mutli-ion plasma more terms 
accounting for the different ions should be included on 
the right hand side of Eq. (36). 

For positive ions, the flux at the electrode is domi- 
nantly due to drift because of the large value of the 
electric field. 



/+ = fj.+n+E" 



(37) 



where E"^ is given by Eq. (11). Setting «+ = at the 
electrode is not appropriate, neither it is desirable from 
the numerical point of view. Although ions are pre- 
sumably completely neutralized at the wall (Auger neu- 
tralization for example), this happens only within a dis- 
tance -10''" m from the wall, which would require 
resolution of a steep ion boundary layer [6]. Negative 
ions are not energetic enough to overcome the potential 
barrier at the wall, hence «_ = 0. For neutrals, a flux 



boundary condition is used. For example for CI atoms 
recombining on the wall with probability yci^cij , one 
has 



7cl-»Cl2 

4 



TTTOcl 



«ci«-7ci+^ci J a*. 



(38) 



where yci+^ci is the probability of ion neutralization on 
the wall, usually assumed unity. This expression states 
that the net flux of CI at the surface equals the flux of 
CI atoms recombining (or in general reacting) on the 
surface minus the flux of CI atoms forming on the 
surface by neutralization of CI* ions. The latter term is 
not important in LDP, but can become very important in 
HDP sources. For highly reactive walls, the concentra- 
tion at the surface may be set equal to zero. For example, 
if metastable atoms are assumed to deactivate with unity 
probability upon striking the wall, «. = 0. Chantry [50] 
has proposed an extrapolation method to specify the 
boundary condition at low pressures. 

The boundary condition for electron mean energy is 
written essentially in the form of an energy balance at 
the electrode [6,19,23]. 



9= = 



kTc 



Sk?t 



Me-H- 



■r+ 



kT^ J, . (39) 



When y+ is zero, combination of Eqs. (36) and (39) with 
the expression for q^ (Eq. ( 1 4)) suggests that the electron 
temperature gradient is zero at the electrode. Some 
workers have assumed a constant electron temperature 
at the electrode (e.g., kT^ = 1 eV). This is convenient 
from the numerical point of view, but the actual value of 
the electron temperature is unknown. In addition, the 
temperature at the wall may be a function of time. 
Kinetic simulations do not have the problem of specify- 
ing a temperature boundary condition; indeed the elec- 
tron energy as a function of time at the wall is an output 
of the simulation [33]. The boundary condition on the 
electric potential is specified as y = on grounded 
walls. On the rf driven electrode the potential is 



V = Vdc + Vrf sin (ot. 



(40) 



The self-bias voltage Vdc has to be found as part of the 
solution. The usual approach is to specify Vrf and then 
adjust the value of Vdc during the simulation so that the 
net charged particle (electrons and ions) current at the 
capacitively coupled electrode is zero [27, 35]. 
Gogolides and Sawin [17] and Dalvie et al. [23] have 
used a current boundary condition on the driven elec- 
trode instead of Eq. (40), assuming singly charged ions. 
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ej+ - e/e + Sa 



dE 
dt ' 



/osin a>t. 



(41) 



Equation (41) implies that the total current, composed of 
the particle (ion and electron) and the displacement cur- 
rents, is forced to be /o sin cot. On insulators the follow- 
ing boundary condition is specified [21, 23], 



eJ+ - e/e + ^0 



dE 
dt 



di dt 



(42) 



which assumes no surface conduction on the insulator. 
Here Vi, di, and ei are the voltage drop across the insula- 
tor, thickness and dielectric constant of the insulator, 
respectively. Eq. (42) is a current continuity equation 
analogous to Eq. (41) and implies that the total current 
(conduction plus displacement current from the plasma) 
equals the displacement current through the insulator. 
Boundary conditions used in fluid simulations have 
been discussed by Wilcoxson and Manousiouthakis [5 1 , 
52]. 

Different sets of initial conditions have been used 
with a varied degree of success. For example, a uniform 
charge density many orders of magnitude lower than the 
expected steady-state value may be specified at f = 0. 
Sometimes parabolic-like charged species distributions 
seem to work better. In any case it is important to satisfy 
the Poisson equation at f = 0; a potential of V = 
everywhere is a common choice. The electron energy 
can be specified as spatially uniform initially. Conver- 
gence is not guaranteed; however, a converged solution 
can be used as an initial condition for another simulation 
at a different set of nearby operating conditions. A for- 
mal continuation scheme can also be applied to conduct 
parametric investigations [53]. Of course, the final peri- 
odic steady-state solution should be independent of the 
choice of initial conditions, unless physically acceptable 
multiple steady states exist. 



5. Computational Aspects 

Glow discharge simulations are computationally in- 
tensive because they are "stiff in both space and time. 
Spatial stiffness results from the fact that rapid changes 
in the dependent variables occur near and inside the 
sheath, which forms naturally over all surfaces in con- 
tact with the discharge. Under high-pressure low-den- 
sity conditions the sheath thickness may be 10-100 
times less than the characteristic discharge dimension, 
depending on pressure, applied voltage, excitation fre- 
quency, etc. The situation is particularly acute in HDP 
reactors which have a sheath thickness of 100-1000 
times less than the discharge dimension. This is because 
high plasma density results in a small (10s of (xm) 



Debye length Aq, and the sheath thickness is only 1 Os of 
Ad. Another problem results from the fact that the 
sheaths are "highly convective" (or "advective") in na- 
ture, meaning that the "drift" current under the influ- 
ence of the electric field [first term on the right hand 
side of Eq. (7)] far surpasses the diffusion current 
[second term on the right hand side of Eq. (7)] . Numer- 
ical handling of highly convective flows in an accurate 
and efficient manner is still an area of active research in 
computational fluid dynamics (CFD) [34]. Tradition- 
ally, upwind differencing has been used in finite differ- 
ence discretizations or the Streamlined Upwind Petrov- 
Galerkin method in finite element approximations 
(SUPG-FEM) [54] . For glow discharge simulations the 
Scharfetter-Gummel exponential scheme, first used in 
modeling of solid state semiconductor devices, has been 
used extensively [55]. Unfortunately, artificial diffusion 
is inevitably introduced by these methods. Flux 
Corrected Transport (FCT) methods [56, 57] are 
designed to minimize artificial diffusion, but specifica- 
tion of the antidiffusive fluxes poses a problem. The 
Donor Cell Method (DCM) [22] and the More Accurate 
FCT (MAFCT) [26] have also been used for plasma 
simulation. 

Stiffness in time is a problem since electron, ion, and 
neutral species response times are <1 ns, 1 (xs to 10 (xs, 
and 10 ms to 100 ms or longer, respectively. The 
smallest operational time scale that needs to be resolved 
is that of the applied excitation frequency, which is typ- 
ically 13.56 MHz (a period of 73.4 ns). Assuming that 
the time step has to be at least 20 times smaller than the 
excitation frequency and that the simulation has to con- 
tinue for 100 ms to capture neutral chemistry, one would 
require some 3X10' time steps to complete the simula- 
tion. Clearly this is a tremendous computational task. To 
make matters even worse, the time step must usually be 
much smaller than that required to resolve the excitation 
frequency. For example to assure stability, explicit time 
integration schemes [58] must have a time step which is 
limited by the Courant-Friedrich-Lewy (CFL) condition 
[36], At<Ax/v^^, where At is the time step size. Ax is 
the grid spacing and i^max is the maximum species veloc- 
ity. Explicit integration is straightforward to execute but 
may require many time steps (> 1 ,000) per cycle to reach 
the periodic steady state. To relax the CFL constraint 
implicit or semi-implicit integration schemes [58] of the 
transport equations (Eqs. (4), (5), and (12)) are usually 
employed [11, 19, 22, 59]. Implicit integration requires 
fewer time steps, but results in a system of complex 
nonlinear equations that may be solved by iterative 
methods. 

For high charge densities, the dielectric relaxation 
(DR) time [22,34,36] imposes severe limitations on time 
step (even more so than the CFL condition) when 
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the Poisson Eq. (23) is integrated explicitly 
AfoR = eo/ie/ji^n^ + e(jb+n+). 



(43) 



For HDP simulations, for example, where the charge 
density is in the range of 10" m"^ to 10'* m~^, one is 
constrained to AfDR<10''^ s! Implicit or semi-implicit 
schemes can overcome this limitation. When Poisson's 
equation, Eq. (23), is solved separately from the charge 
continuity equations (either implicitly or explicitly) the 
time adwncement of the potential is explicit in nature as 
the charge density is evaluated at time t. Therefore, this 
approach is bound to advance in time with at most A^dr- 
However, Ventzek et al. [22] devised a semi-implicit 
update technique for the electric potential that allowed 
the time steps to exceed Af^R by 1-2 orders of magni- 
tude. The method is based on approximating the charge 
density at a future time, p(t + At), by relying on the 
present values of p(t) and its derivative (dp/dt), (i.e., 
p(t + At) ~ p{t) + dpldt),At). Stewart et al. [59] over- 
came the dielectric relaxation time step limitation by 
removing the strong coupling between the electron con- 
tinuity and Poisson's equation. This was achieved by 
substituting the Laplacian of the electric potential (i.e., 
V^y) appearing in the electron continuity equation (Eq. 
(4) written for electrons) by the right hand side of 
Eq. (23). 

Formal acceleration schemes based on the Newton- 
Raphson method [17, 19, 53] or heuristics based on 
extrapolation [22, 60-62] have been used to accelerate 
convergence to the periodic steady state. It has been 
estimated that these techniques can speed up conver- 
gence by many orders of magnitude. 

In order to decouple the disparate time scales of elec- 
tron and neutral transport, Kushner and co-workers [22, 
62] and Lymberopoulos and Economou [19, 20] used a 
modular approach which can be thought of as an equa- 
tion splitting technique. An example, used for a capac- 
itively-coupled argon discharge, is shown in Fig. 2. The 
glow discharge module includes the electron and ion 
density continuity [Eq. (4)] and the drift-diffusion [Eq. 
(7)] equations which are solved simultaneously with the 
Poisson equation [Eq. (23)]. The electron energy (tem- 
perature) equation [Eq. (13)] is then updated using the 
new charge densities and fields. In turn, the neutral 
metastable density equation [Eq. (4)] is solved on a 
staggered mesh which does not have to be as fine as the 
grid used for the glow discharge and electron tempera- 
ture modules. The simulation then returns to the glow 
discharge module with updated values of the electron 
energy and metastable density. The equations for the 
effective electric field to which the ions respond [Eq. 
(11)] are also solved periodically. Information is thus 
cycled back-and-forth among the modules until conver- 



gence. This approach can be extended to the GEC-ICP 
cell by including a module that solves for the azimuthal 
electric field [Eq. (31)]. This has been done for simulat- 
ing the GEC cell [62] and other HDP reactors [46]. 
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Fig. 2. Modular approach used to perform a 2-D plasma simulation 
of an argon rf discharge. 



6. Results and Discussion 

6.1 GEC/CCP Cell 

Although there have been quite a few 1-D simulations 
of the GEC-CCP cell (see paper by Govindan and 
Meyyappan in this Special Issue), the authors are aware 
of only three published 2-D simulations of the GEC/ 
CCP cell [20, 21, 27], and one published simulation [62] 
of the GEC-ICP cell. Preliminary results of a hybrid 
2-D simulation of the GEC-CCP cell were reported by 
Pak and Riley [63]. Other works in two-dimensional 
geometries that are relevant to the GEC-CCP cell are the 
papers of Dalvie et al. [23], and Wu and co-workers 
[26]. Still more 2-D capacitively-coupled rf glow dis- 
charge simulations have been reported by Tsai and Wu 
[24], Paschier and Goedheer [35], and Wilcoxson and 
Manousiouthakis [52], all using the fluid approximation, 
and by Vahedi et al. [30] who used a PIC-MCC 
approach. 

Lymberopoulos and Economou [20, 21] simulated an 
argon discharge including the transport of neutral 
metastables in a self-consistent manner. The problem 
was simplified by considering only one (lumped) 
metastable state, and the relatively simple chemistry 
shown in Table 1 . They used a 0-D Monte Carlo simula- 
tion to express the electron-impact reaction rate coeffi- 
cients as a function of mean electron energy [19]. One 
of the results of that work was the significance of 
metastable neutral atoms in sustaining the discharge by 
a two-step ionization process (i.e., excitation to the 
metastable level first followed by ionization of that level. 
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see Reactions 1 and 3 in Table 1). The importance of 
two-step ionization was recognized earlier in connection 
with dc discharges (positive column) [64, 65]. Figure 3 
shows a comparison of the time-average electron den- 
sity obtained by a 1-D rf (13. 56 MHz) simulation of an 
argon discharge with (solid line) and without (dashed 
line) metastables [19]. Including metastables results in 
several times higher electron density which is then in 
accord with experimental measurements [66]. Metasta- 
bles have an ionization threshold of only 4.14 eV as 
compared to 15.7 eV for direct ionization from the 
ground state. Although the metastable density is orders 
of magnitude smaller than the ground state species, the 
corresponding ionization coefficient is orders of magni- 
tude greater; in fact it is such that the two-step ionization 
is the dominant mechanism for electron production 
under these conditions (133.3 Pa pressure). The impor- 
tance of metastables in sustaining the discharge points to 
the possible effects that minute amounts of impurities 
(e.g., air or moisture from a vacuum leak) can have, 
since these impurities can quench metastables very 
effectively. 




T 

0.2 0.4 0.6 0.8 1 

Dimensionless Position, ^, 

Fig. 3. Electron density between the electrodes of a capacitively 
coupled rf argon discharge. Solid line: including metastables in the 
calculation in a self-consistent manner. Dashed line: without metasta- 
bles. Conditions: 133.3 Pa, 300 K, 100 V peak-to-peak, 13.56 MHz 
rf voltage (from Ref. [19]). 



The time-average electron density distribution from a 
2-D rf (13.56 MHz) simulation of an argon discharge at 
a pressure of 133.3 Pa is shown in Fig. 4 [20]. A push- 
pull scheme was used to power the Cell which resulted 
in a symmetric discharge (no dc bias). Only part of the 
Cell of Fig. la is shown for clarity. Near the reactor 
centerline the electron density has a cosine-like profile 



similar to that of the 1-D simulation of Fig. 3. A peak in 
electron density is observed in the radial direction. The 
electron density drops off drastically beyond the elec- 
trode edge, i.e., the plasma is well confined. This is 
because of the relatively high pressure (133.3 Pa) low 
voltage operation (60 V peak-to-peak) and the fact that 
the discharge is symmetric. 

The radial peak in charge density has been attributed 
to the synergistic effect between the axial and radial 
electric fields near the edge of the electrodes. A radial 
(space charge) electric field develops as the more mo- 
bile electrons escape the plasma between the two elec- 
trodes into the surrounding chamber. The radial electric 
field extends over a region ~ 1 cm inwards (towards the 
centerline) from the electrode edges (Fig. 5). Away from 
the edge the electric field is directed axially and heats 
the electrons primarily near the plasma-sheath interface. 
Near the edge, the electrons acquire extra energy due to 
the radial electric field and hence enhance the ionization 
resulting in a radial peak in the charge density. Radial 
peaks in electron and ion density have been measured in 
the GEC-CCP Cell by Overzet and Hopkins [67]. 

The enhanced electron energy near the electrode 
edges is also reflected in the neutral metastable density 
profiles: hot spots in metastable density develop near the 
edges as seen in Fig. 6. The metastable density is de- 
pressed along the midplane between the electrodes be- 
cause the main loss mechanism for metastables is 
quenching by electrons (reaction 5 in Table 1); and the 
electron density peaks along the midplane. A structure 
similar to that shown in Fig. 6 has been observed exper- 
imentally by Greenberg and Hebner [66] in the Cell for 
a corresponding helium metastable state (see Fig. 7). 

The importance of metastables in contributing to ion- 
ization at a pressure of 133.3 Pa was discussed in con- 
nection with Fig. 3. Figure 8 shows the electron produc- 
tion rate due to direct ionization (Reaction 2 in Table 1) 
and two-step ionization (reaction 3 in Table 1), for a 
pressure of 13.3 Pa. Two-step ionization is only 10 % of 
the direct ionization in this case. As the pressure is 
lowered, the mean electron energy increases. Higher 
electron energies favor more endothermic reactions 
(direct vs two-step ionization). Even at this lower pres- 
sure metastables make a significant contribution to ion- 
ization. Ventzek et al. [62] have shown that the two-step 
ionization continues to be significant even down to 1.33 
Pa. Apparently even at 1.33 Pa in the GEC-ICP 
discharge [62], the electron energy is not high enough 
for direct ionization to completely dominate. Another 
reason is that the metastable density Mat* depends rela- 
tively weakly on pressure; it appears to change by only 
a few times as the pressure changes by a factor of 100 
(compare results of Ref. [19], [21], and [62] for pres- 
sures of 133.3 Pa, 13.3 Pa, and 1.33 Pa respectively). 
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Fig. 4. Electron density (in 10 m ) in an argon discharge in the GEC Cell of Fig. 
la for a symmetric (push-pull) rf drive. Conditions: 133.3 Pa, 300 K, 60 V peak-to- 
peak, 13.56 MHz rf voltage (from Ref. [20]). 
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Fig. 5. Radial electric field distribution in the GEC Cell of Fig. la for symmetric (push-pull) 
rf drive. Conditions: 133.3 Pa, 300 K, 60 V peak-to-peak, 13.56 MHz rf voltage. 
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Fig. 6. Argon metastable density in the GEC Cell of Fig. la for symmetric (push-pull) rf drive. 
Conditions: 133.3 Pa, 300 K, 60 V peak-to-peak, 13.56 MHz rf voltage. 
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Fig. 7. Measured helium metastable (2 'X) density in the GEC Cell. 
Conditions: 133.3 Pa, 300 V peak-to-peak, 13.56 MHz rf voltage 
(from Ref. [66]). 



Thus, the ratio of metastable to ground state density 
increases with decreasing pressure and that makes the 
two-step ionization contribution significant. 

Boeuf and Pitchford [27] used an approach very sim- 
ilar to that of Lymberopoulos and Economou [20, 21] 
except that Boeuf and Pitchford did not include 
metastable neutral transport. Their method of solution 
was based on a finite difference approximation of the 
spatial derivatives using the Scharfetter-Gummel expo- 
nential scheme. They used a uniform 41X41 grid. For 
an electrode spacing of 2.54 cm, this corresponds to a 
grid spacing of 0.6 mm between the electrodes. An 
explicit time integration was used to advance the simula- 
tion to the periodic steady-state reached after some 1000 
rf cycles (13.56 MHz) with 500 time steps per cycle. 
They compared their simulation results with the data of 



Overzet and Hopkins [67] and they found good agree- 
ment with those measurements (Fig. 9). The agreement 
was much better at 33.3 Pa (shown in Fig. 9) than at 13.3 
Pa which is surprising since Boeuf and Pitchford did not 
include metastables in their calculation and the 
metastable contribution is expected to be more signifi- 
cant at the higher pressure. 

The time-average radial electron density profiles at a 
distance of 1.25 cm from the powered electrode are 
shown in Fig. 10a as a function of applied rf voltage 
amplitude, and in Fig. 10b as a function of pressure [27]. 
It is clear that higher voltages and/or pressures enhance 
the off-axis radial peak in electron density. As voltage 
increases the radial electric field increases as well which 
in turn increases the ionization rate near the radial edge 
of the plasma. Also, as the pressure is lowered diffusion 
becomes more facile smoothing the concentration 
gradients. 

"ibung and Wu [26] simulated a 13.56 MHz helium 
discharge in a geometry relevant to the GEC-CCP cell. 
They truncated the surrounding chamber by placing a 
cylindrical solid wall confining the discharge to a radius 
of 5.08 cm (2 in). The electrode spacing in their simula- 
tion was 2.54 cm (1 in) as in the Reference Cell. They 
used a fluid approximation with the full momentum and 
energy equations (Eqs. (5) and (12)) for the electrons 
and the drift-diffusion approximation with an effective 
electric field (Eqs. (7) and (11)) for ions. They did not 
consider metastable transport in their simulation. 
Electron-impact reaction rate coefficients and transport 
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Fig. 8. (a) Rate of ionization of argon by the direct mechanism (reac- 
tion 2 of Table 1). Contour values should be multiplied by lO" m"^ s"'. 
(b) Rate of ionization of argon by the two-step mechanism (reaction 
3 of Table 1). Contour values should be multiplied by lO'* m"^ s"'. 
Conditions: 13.3 Pa, 300 K, 60 V peak-to-peak, 13.56 MHz rf 
voltage. 



properties were obtained by a 0-D DC Monte Carlo 
simulation. Figure 11 shows that a peak in the radial 
profile of the time-average ion density appears in this 
case as well. Young and Wu attributed this peak to the 
presence of a radial electric field which provides extra 
heating of the electrons in the region close to the radial 
wall. In this case the radial electric field is present be- 
cause of the sheath which forms naturally over the radial 
wall. It is interesting to note that their result (i.e., off- 
axis radial peaks in charge density) is similar to that 
obtained in the open GEC-CCP geometry of Fig. 1 
where the radial confining wall is away from the elec- 
trodes [20,21,27]. 

Dalvie et al. [23], also observed an off-axis radial 
peak in the charge density in a 2-D simulation of a 
reactor with a radial confining wall. They used the fluid 
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Fig. 11. Time-average ion density in a helium 13.56 MHz rf discharge at 133.3 Pa and 
150 Vcm"' (vacuum) field between the electrodes. 



equations with the drift-diffusion approximation (no 
effective field for ions), and Arrhenius rate expressions 
for the ionization and excitation rates of argon as a 
function of electron temperature. The authors did not 
consider metastable transport. They used a sinusoidal 
total current boundary condition at the powered 
electrode (Eq. (41) with /o = 0.175 A, and co = 2 ttv, 
v= 13.56 MHz). The radial wall was an insulator 
(quartz) which was grounded on the opposite side of 
that exposed to the plasma. When the insulator was 
"thick," the discharge was nearly symmetric since the 
resistance was too high for appreciable (displacement) 
current to flow through the insulator. The time-average 
ionization rate at a pressure of 66.7 Pa is shown in Fig. 
12. Local maxima in ionization are seen near the dis- 
charge corners which the authors attributed to extra 
electron heating at the corners due to the focusing of the 
rf current by the radial sheath. The focusing of the 
current is seen in Fig. 13 which displays the root-mean- 
square current density in the cell. As in the case of the 
simulations of Young and Wu [26] and Lymberopoulos 
and Economou [20, 21], the presence of a radial sheath 
contributes to extra electron heating. Since the spatial 
distribution of ionization and excitation rates in argon 
are similar (both are high threshold processes), the local 
maxima in Fig. 12 would produce a local enhancement 
in the density of metastable species, should the authors 
have included metastables in their calculation. This 
would be consistent with Fig. 6. 



Anode 







Cathode (cm) 

Fig. 12. Time-average ionization rate in an argon discharge at 66.7 
Pa and 0.175 A, 13.56 MHz rf peak current (see Eq. (41)). Contour 
values increase in steps of 2X10^ m"^' s"' (from Ref. [23]). 



6.2 GEC/ICP Cell 

Kushner and coworkers [62] have developed a com- 
prehensive 2-D simulation of the GEC-ICP cell shown 
in Fig. 1(b). They used a modular approach which is an 
extension of their previous work on 1 -D glow discharge 
simulations [18] and which is in the same spirit of the 
modular approach (Fig. 2) followed by Lymberopoulos 
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Fig. 13. Root-mean-square current density in an argon discharge at 66.7 Pa and 0.175 A, 13.56 MHz rf peak current 
(see Eq. (41)) (from Ref. [23]). 



and Economou [20] . Fluid equations were used to com- 
pute the electron, ion and neutral species densities. An 
electromagnetics module solved for the azimuthal elec- 
tric field distribution [Eq. (31)] and Eq. (18) was used 
to obtain the magnetic fields. Poisson's equation was 
solved implicitly simultaneously with the charged spe- 
cies density continuity equations. The electromagnetic 
fields were used in a 2-D Monte Carlo simulation per- 
formed in regular intervals during the progress of the 
overall simulation to determine the time-average elec- 
tron transport properties and rate coefficients of elec- 
tron-impact reactions. Also, a Navier-Stokes hydrody- 
namics module was used to calculate the gas velocity 
distribution. Information was cycled among the modules 
until convergence. This is a hybrid simulation in the 
sense that the electron impact rate coefficients are cal- 
culated by a kinetic Monte Carlo scheme and not by a 
(fluid) energy equation. 

The results of the GEC-ICP simulation by Ventzek et 
al. [62] are shown in Figs. 14 and 15. Figure 14 shows 
the time-average electron density (right) and plasma 
potential (left). The stove-top coil dimensions are also 
shown in Fig. 14. The plasma is rather well confined 
despite the low operating pressure of 1. 33 Pa. This is 
consistent with experimental measurements in the GEC- 
ICP cell. The electron density drops by an order of 



magnitude a small distance beyond the edge of the coil, 
and decays to very small values in the surrounding 
chamber. The plasma potential peaks at around 20 V. 
The negatively charged alumina support at the lower 
part of the reactor is also shown. The electrons are fairly 
hot under the coils (Fig. 14, left) but they cool off 
drastically in the surrounding chamber. The ionization is 
confined in the main plasma volume under the coils 
(Fig. 15) since both the electron density and tempera- 
ture drop off sharply as a function of radius beyond the 
coil edge. Good agreement was obtained for the electron 
density as a function of power between the computed 
and measured values [68]. 

Economou and B artel [69] have used a kinetic Direct 
Simulation Monte Carlo (DSMC) method to simulate 
gas flow and pressure distribution in the GEC-ICP cell. 
DSMC is expected to be more accurate than a fluid 
simulation at pressures down to 0. 13 Pa at which the 
mean free path of species exceeds the reactor spacing 
(Kn > 1). Figure 16 shows the distribution of the radial 
component of the convective flow velocity (Fig. 16a) 
and the corresponding pressure distribution (Fig. 1 6b) in 
only part of the cell of Fig. lb. The flow rate was 140 
seem and plasma effects were not included in this 
simulation (neutral flow only). The inflow port was 
modeled as a point source with choked flow at the inlet. 
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Fig. 14. Time-average electron density (right) and potential (left). 1.33 Pa inductively-coupled argon plasma in 
the GEC Reference Cell of Fig. lb (from Ref. [62]). A contour value of 100 corresponds to 4X10" m"^' electron 
density. 
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Fig. 15. Time-average rate of ionization (right) and electron temperature (left). 1.33 Pa inductively-coupled 
argon plasma in the GEC Reference Cell of Fig. lb (from Ref. [62]). 



490 



Volume 100, Number 4, July- August 1995 

Journal of Research of the National Institute of Standards and Technology 



Radial Velocity (m/s) 

-213 -169 -125 -82 -38 6 



^ 0.04 



2 




0.00 
0.00 



0.04 0.06 0.08 

RADIAL POSITION 



(m) 

0.06 



0.00 
0.00 



Pressure (Pa) 




0.02 



0.04 0.06 

RADIAL POSITION 



0.08 



0.10 



(m) 
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One observes that the average flow velocity can exceed 
250 m s"' under these conditions. The "plume" of the 
inflow jet is clearly seen. Also, the flow is relatively very 
slow in the central region of the Cell. There are sub- 
stantial pressure gradients in the system. A region of 
~ 5X10'^ m radius around the inlet port has pressures 
exceeding 1.33 Pa. The main body of the cell is at 
around 0.4 Pa, and the exhaust is at less than 0.13 Pa. 
Detailed DSMC simulations of plasma flow [3 1 ] have 
shown that substantial total density gradients can exist in 
the reactor due to gas heating (by charge exchange and 
the Frank-Condon effect) and ion pumping. Ion pump- 
ing refers to the phenomenon in which neutrals are 
ionized in the bulk of the plasma and the resulting ions 
are driven to the reactor walls by the space charge fields 
where they recombine turning back into neutrals. Thus 
the walls are populated with more neutrals compared to 
the bulk plasma. Plasma-DSMC calculations of the 
GEC-ICP Cell are currently in progress. 



7. Summary and Outlook 

During the past few years, multidimensional self -con- 
sistent plasma simulations which account for the cou- 
pled effects of charged and neutral species transport and 
chemistry have been developed. These simulations vary 
in their degree of detail from kinetic to fluid to hybrid 
simulations. Also, different degrees of approximation 
are used within the same group (e.g., fluid) of simula- 
tions. Detailed comparisons with experimental data are 
necessary to decide which degree of approximation is 
adequate for accurate determination of some important 
quantities such as the species density profiles, and the 
radical and ion flux and energy uniformity along the 
electrodes. For example, how low in pressure can the 
fluid approximation be used, and what is the range of 
pressure and frequency for which the drift-diffusion 
approximation is correct? Detailed comparisons with 
data is now actively pursued by many research groups. 

For complex chemical systems the accuracy of the 
simulation may be limited by the lack of knowledge of 
cross sections for electron-impact reactions and plasma 
chemistry. In etch or deposition plasmas, knowledge of 
surface chemistry is viewed as an even more important 
limitation. This becomes more acute as the operating 
pressure decreases. Experiments in well characterized 
systems such as the GEC reference cell, in combination 
with plasma simulation, will continue to enhance our 
fundamental knowledge about the plasma dynamics. 
Many more comparisons with experimental data are 
needed to "tune" the models and provide simulation 
tools with predictive capabilities. Considering that 2-D 
self-consistent rf plasma simulations are only a few 



years old, such activity is expected to be vigorous in the 
near future. 

The near future will also witness further use of the 
GEC-ICP cell as high density plasma operation be- 
comes even more important industrially. Also, more 
studies will be conducted with reactive gas chemistries 
for wafer etching. The development of diagnostics (es- 
pecially non-intrusive optical diagnostics) and sensors 
for real time process control will continue. Although at 
an early stage of development multidimensional plasma 
simulations have shown great promise in reproducing 
many of the dominant features observed experimentally. 
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